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We evaluate the non-Markovian finite-temperature two-time correlation functions (CF's) of system 
operators of a pure-dephasing spin-boson model in two different ways, one by the direct exact 
operator technique and the other by the recently derived evolution equations, valid to second order 
in the system-environment interaction Hamiltonian. This pure-dephasing spin-boson model that is 
exactly solvable has been extensively studied as a simple decoherence model. However, its exact 
non-Markovian finite-temperature two-time system operator CF's, to our knowledge, have not been 
presented in the literature. This may be mainly due to the fact, illustrated in this article, that 
in contrast to the Markovian case, the time evolution of the reduced density matrix of the system 
(or the reduced quantum master equation) alone is not sufficient to calculate the two-time system 
operator CF's of non-Markovian open systems. The two-time CF's obtained using the recently 
derived evolution equations in the weak system-environment coupling case for this non-Markovian 
pure-dephasing model happen to be the same as those obtained from the exact evaluation. However, 
these results significantly differ from the non-Markovian two-time CF's obtained by wrongly directly 
applying the quantum regression theorem (QRT), a useful procedure to calculate the two-time CF's 
for weak-coupling Markovian open systems. This demonstrates clearly that the recently derived 
evolution equations generalize correctly the QRT to non-Markovian finite-temperature cases, ft is 
believed that these evolution equations will have applications in many different branches of physics. 
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I. INTRODUCTION 

A quantum system is inevitably subject to the influ- 
ence of its surroundings or environments pHa] ■ An envi- 
ronment usually consists of a practically infinite number 
of degrees of freedom and acts statistically as a whole 
identity referred as a reservoir or bath of the open quan- 
tum system. Most often, one is concerned with only the 
system dynamics and the key quantity is the reduced sys- 
tem density matrix p(t) defined as the partial trace of the 
total system-plus- reservoir density operator px(t) over 
the reservoir degrees of freedom; i.e., p(t) — Tr r[pt(£)]- 
If the time evolution of the reduced density matrix that 
can be Markovian or non-Markovian is known, one is able 
to calculate the (one-time) expectation values or quan- 
tum average of the physical quantities of the system op- 
erators. But knowing the time evolution of the reduced 
density matrix is not sufficient to calculate the two-time 
(multiple-time) correlation functions (CF's) of the sys- 
tem operators in the non-Markovian case 0-Q . 

In the Markovian case, an extremely useful proce- 
dure to calculate the two-time (multiple-time) CF's is 
the so-called quantum regression theorem (QRT) [l|-|4| 
that gives a direct relation between the time evolution 
equation of the single-time expectation values and that 
of their corresponding two-time (multiple-time) CF's. So 
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knowing the time evolution of the system reduced den- 
sity matrix allows one to calculate all of the two-time 
(multiple-time) Markovian CF's. For the non-Markovian 
case, it is known that the QRT is not valid in general [lOl — 
fl3j |. Recently, using the stochastic Schrodinger equation 
approach and the Heisenberg equation of system opera- 
tor method, an evolution equation, valid to second order 
in system-environment coupling strength, for the two- 
time (multiple-time) CF's of the system operators has 
been derived for an environment at the zero temperature 
and for a system in an initial pure state [7H9|. This evo- 
lution equation has been applied to calculate the emis- 
sion spectra of a two-level atom placed in a structured 
non-Markovian environment (electromagnetic fields in a 
photonic band-gap material) [14J. In Ref. [8j], an evolu- 
tion equation for the reduced propagator of the system 
state vector, conditioned on an initial state of the en- 
vironment differing from the vacuum, was derived using 
the stochastic Schrodinger equation approach. It is thus 
possible to use the reduced propagator to evaluate the 
expectation values and CF's of the system observables 
for general environmental initial conditions, not neces- 
sarily an initial vacuum state for the environment Q. 
By using another commonly used open quantum system 
technique, the quantum master equation approach [lH6j , 
we are able to extend the two-time CF evolution equa- 
tion to a non-Markovian finite-temperature environment 
for any initial system-environment separable state. The 
detailed derivation will be presented elsewhere [15| but 
the essential results will be summarized in Sec. |TTJ The 



derived evolution equation that generalizes the QRT to 
the non-Markovian finite-temperature case is believed to 
have applications in many different branches of physics. 

The purpose of this article is twofold: (a) We show 
that in general the time evolution of the reduced den- 
sity matrix of the system (or the reduced quantum mas- 
ter equation) alone is not sufficient to calculate the two- 
time CF's of the system operators of non-Markovian open 
systems, even in the weak system-environment coupling 
case. We present an evaluation of an exactly solvable 
non-Markovian model, i.e., a pure-dephasing spin-boson 
model [H, Il6l - l2l| . to justify the statement. The exact 
non-Markovian finite-temperature two-time CF's of the 
system operators of this model, to our knowledge, have 
not been presented in the literature, (b) This exactly 
solvable model allows us to test the validity of the derived 
non-Markovian finite-temperature evolution equation of 
two-time CF's presented in Sec. [TTJ It will be shown that 
the two-time CF's obtained using the evolution equa- 
tion in the weak system-environment coupling limit [la ] 
in Sec. [TT| for the exactly solvable non-Markovian model 
happen to be the same as those obtained from the ex- 
act evaluation. However, these results significantly differ 
from the non-Markovian CF's obtained by wrongly ap- 
plying directly the QRT. This demonstrates clearly that 
the derived evolution equations generalize correctly the 
QRT to non-Markovian finite-temperature cases. 

The article is organized as follows. We first summa- 
rize the important results of the newly obtained evolu- 
tion equations [15| that generalizes the QRT to the non- 
Markovian finite-temperature case in Sec. [TTJ After brief 
description of the pure-dephasing spin-boson model in 
the beginning of Sec. lIIIl we calculate the exact time evo- 
lution of the reduced density matrix of the system and 
one-time expectation values in Sec. lIII Al The exact two- 
time CF's are evaluated in subsection IIIIBI In Sec. IIV1 
we use the derived evolution equations in Ref. [15| to cal- 
culate the one-time and two-time CF's. It is shown that 
the results obtained in Sec. IIVI are the same as those 
by the exact evaluation in Sec. IIIII This demonstrates 
the validity and practical usage of the derived evolution 
equations in Ref. [l^- Numerical results and discussions 
are presented in Sec. [V] A short conclusion is given in 
Sec. ED 



II. EVOLUTION EQUATION OF 

NON-MARKOVIAN FINITE-TEMPERATURE 

TWO-TIME CF'S 

A class of systems considered in pHS! is modeled by 
the Hamiltonian 

H = Hs + Hj + H r 

= H S + J2 h 9\ (tia x + Laty + Y, ^4 a A> (1) 

A A 

where Hs and Hr are system and environment Hamil- 
tonians, respectively, and Hj stands for the Hamiltonian 



that describes the interaction between the system and 
the environment. So L acts on the Hilbert space of the 
system, a\ and a\ are creation and annihilation opera- 
tors on the environment Hilbert space, and g\ and u\ 
are the coupling strength and the frequency of the Ath 
environment oscillator, respectively. The derivations of 
the non-Markovian evolution equations of the two-time 
(multitime )CF's for the general Hamiltonian model ([1]) 
in Refs. 041| (Eq. (6) in Ref. Q, Eq. (31) in Ref. @ and 
Eq. (60) in Ref. Q) are presented for an environment at 
the zero temperature and for a system state in an initial 
pure state. It was mentioned in Ref. Q that it is possible 
to use the reduced stochastic system propagator that cor- 
responds to an initial state of the environment different 
from the vacuum to evaluate the single-time expectation 
values and multitime CF's with more general initial con- 
ditions. But only a master equation that is conditioned 
on initial bath states and is capable of evaluating the 
single-time expectation values of system observables for 
general initial conditions, both for an initial pure state 
and mixed state, was derived @. In Refs. [iHil, calcu- 
lations of the two-time CF's of system observables for 
dissipative spin-boson models in thermal baths are, how- 
ever, presented even though in their derivations of the 
two-time (multitime) evolution equations, the bath CF's 
are given in its zero-temperature form. This is possible 
due to the reason that for a system-environment model 
with a Hermitian system operator L = L^ coupled to 
the environment, the linear finite-temperature stochas- 
tic Schrodinger equation could be written in a simple 
form of the zero-temperature equation [19L |22I | if the zero- 
temperature bath CF is replaced with its corresponding 
effective finite-temperature bath CF. As a result, the evo- 
lution equation of thermal two-time (multitime) CF's for 
a Hermitian coupling operator L = L^ also becomes equal 
to its zero-temperature counterpart with the replacement 
of the zero-temperature bath CF with its effective finite- 
temperature bath correlation kernel. It is for this reason 
that the dissipative spin-boson model with a thermal en- 
vironment can be studied with the two-time (multitime) 
evolution equations derived in Refs. [7H9J], since in that 
model L = a x = L^ . But this reduction of the finite- 
temperature evolution equation to its zero-temperature 
form 7 9] is not valid for more general non-Markovian 
finite-temperature cases where the system coupling op- 
erators are not Hermitian, i.e., L ^ L* . In other words, 
if the system operator coupled to the environment is not 
Hermitian L ^ L\ the two-time (multitime) differential 
evolution equations presented in Refs. [IH9| are valid for 
a zero-temperature environment only. 

By using another commonly used open quantum sys- 
tem technique, the quantum master equation approach 
QHa], it is possible to obtain in the weak system- 
environment coupling limit a two-time evolution equa- 
tion for non-Markovian finite-temperature environments 
with both Hermitian and non-Hermitian system coupling 
operators and for any initial system- environment sepa- 
rable states. The detailed derivation will be presented 



elsewhere [151 ] but the important results are summarized 
here. The second-order evolution equations of the single- 
time expectation values for the class of systems modeled 
by the Hamiltonian (JTJ) is 

diAitx^/dh 
= (i/h)Tr s ({[H s ,A}}(h)p(0)) 

drTrs 



(a*(ti-r){lt( r _t0[4£]}(*i)p(0) 

+a(h - t) {[L\A]L(t - ti)} (t!)p(0) 
+(3*(t 1 -T){L(T-t 1 )[A,L^}}(t 1 )p(0) 
+ p(t 1 -T)[[L,A}D(r-t 1 )}(t 1 )p(0) 

and that of the two-time CF's can be obtained as 
d(A(h)B(t 2 ))/dti 
= (i/h)Tr s ({[H s ,A]}(h)B(t 2 )p(0)) 

+ f ' drTr s 
Jo 

(a*(h - t) {D(t - ti)[A, L]j (ti)B(*a)p(0) 
+o(ti - r) [\L\A]L(t - ti)} (*i)B(ta)p(0) 
+/?*(*! - r) {X(r - t!)[A, Lt]} (i 1 )B(t 2 )p(0) 
+ /3(t! - r) {[L, A]£t( T - tl )} (*i)B(«a)p(0) 



(2) 



drTrs 

^-^{[it.AlXtijIlB.ICr-taJlJ^MO) 
/?(*! - r) {[£, A]} ( tl ) {[£, £t(r - fc)]} (t a )p(0)) . (3) 



Here L(t) = exp (iHgt/h) Lexp (—iHst/h) is the system 
operator in the interaction picture with respect to H$, 
and 

«(r - s) = £> A + l)l3A| 2 e'^ (T - s) , (4) 

A 

£(t-s) = 5> a |<7a|V^ t - s >. (5) 

A 

are known as the environment CF's: <x(t — s) = 

EaS^aC"^^'^'"!'^)) and P( T ~ s ) 

T,\9\a{(T)J2x 9\>a\>(s)}, where 6a(t) = a x er lu >- T 

and OaC 7 ") = a \ elu)XT are the reservoir operators in the 
interaction picture. 

We note here that for a Hermitian coupling operator 
L = iJ the finite-temperature evolution equations © 
and (|3J) reduce, respectively, to their zero-temperature 
counterparts but with the effective bath CF given by 
a(ti — r) + j3(ti — t) 7-9]. This was pointed out to 
occur in general for iV-time CF's in Refs. pHll- 



III. EXACT EVALUATIONS OF PURE 
DEPHASING SPIN-BOSON MODEL 



Here we consider an exactly solvable pure dephasing 
model of 



Ha 



(hujs/2)a z , 



L = a, 



V 



(6) 



to test the evolution equations © and ([3]). This pure 
dephasing spin-boson model in which [Hs,L] = has 
been extensively studied as a simple decoherence model 
in the literature [8|, [l6l - [21j . But most of the studies focus 
on the discussion of the time evolution of the reduced 
density matrix of the spin, or other one-time expecta- 
tion values of the spin system operators. Recently, the 
two-time CF's of the system operators at the zero tem- 
perature for this model was reported in Ref. [8}. Never- 
theless, to demonstrate the validity and practical usage 
of the finite-temperature non-Markovian evolution equa- 
tion of the two-time CF's ([3]), we present a detailed eval- 
uation of the exact finite-temperature two-time CF's for 
this simple model. These exact non-Markovian finite- 
temperature two-time CF's of the system operators, to 
our knowledge, have not been presented in the literature. 



Reduced density matrix and one-time 
expectation values 



Before we derive the two-time CF's, we evaluate the 
exact time evolution of the reduced density matrix and 
one-time expectation values for the non-Markovian spin- 
boson model. In the interaction picture, the total density 
matrix of the combined (spin plus bath) system at time 
t is given by 



Pt 



(t) = U(t)p T (0)U\t), 



(7) 



where the time evolution operator is 



U(t) 



_ e iH t/h e -iHt/h 



(8) 



Here H = H s + H R , Hj (t) 

exp (iHot/K) Hjexp (— iHgt/h) and T is the time- 
ordering operator which arranges the operators with the 
earliest times to the right. ^From Eqs. ([T]) and ([5]), a 
simple calculation gives 



H I (t) = J2 n 9\<?z( 



e iWA *a A 



e-*"**OA 



(9) 



This result allows us to calculate the time evolution op- 
erator to be (see Appendix A for details) 



result, we obtain 



U(t) = cxp 



Jo A v 

cxp U J\r jy s J2\9x\ 2 e^ {T - s) ) 

exp [ — / dr I ds y 

\ Jo Jo A 



gx f e - iux(T -s) (1Q) 



The time integrations in the exponents in Eq. (| 10|) can 
be easily and analytically carried out. But we keep them 
in those forms in Eq. (fTU|) so it will be easier to identify 
them with the results in Ref. [15|]. If the time-ordering 
operation in Eq. ([8]) for U(t) were not performed, one 
could have just obtained the first term (line) of Eq. (fTU)) 
for U{t). Thus the second and third terms (lines) of 
Eq. (fTU|) can be considered as the correction terms due 
to the time-ordering operation. 

The reduced density matrix can be obtained 
by tracing over the reservoir's degrees of freedom: 
p(t) = Ti ji[pT(t)]. Suppose initially the state 
Pt(0) — P~t(0) = p(0) <& Ro is factorized, where 
p(0) and Rq are initial system and thermal reser- 
voir(environment) density operators, respectively, and 
R = exp{-H R /k B T)/Tr R [exp(-H R /k B T)]. Then the 
reduced density matrix elements in the interaction pic- 
ture can be written as 

~Pmn(t) = /w(0)Tr fi (u^{t )U {m} (t)R ) , (11) 

where p mn {t) = (m\p{t)\n), U^(t) = (n\U(t)\n), 
in, n — 0, 1 and the states of the two-level system are de- 
fined as er 2 |0) = |0), <r z |l) = -|1). To evaluate Eq. ([IT]), 
the well known formula of 



e A e B =e A+B e^ A ' B K 



(12) 



valid for operators A and B both commuting with the 
commutator [A, B], can be used to combine the evolution 
operators together. One then obtains 



Jo , v 



e~^ T a x 



.(13) 



Then a useful identity [23] for the average over the ther- 
mal reservoir (environment) density operator, Ro, can be 
employed: 



,Ea c\a\+d\a\ 



e kT,xCxd > ,(2n x + l) 



(14) 



where c\ , d\ are complex numbers, and fi\ = 
[exp(Tiuj\/kBT) — 1] _1 stands for the thermal mean oc- 
cupation number of the environment oscillators. As a 



Trr 



U^ 0} (t)U {1} {t)R 



Tr 



UHV{t)U {0} (t)Ro 
ex p(~/ drD(T)\ (15) 



where 



D(t) = 2 / ds[a cB (r - a) + a* cS {r - s)}, (16) 

Jo 



a e s{r - s) = a(r - s) + /3(t - s), 



(17) 



and a(t — r) and (3(t — r) are defined in Eqs. (g]) and (J5JI, 
respectively. It is easy to show that U^ n ^ (t)U^ (t) — 

1. Thus, using 



J and Tr 



R 



UH n }(t)UW(t)Ro 

these results for the reduced density matrix elements 
Eq. (|lip in the interaction picture and then trans- 
forming them back to the Schrodinger picture pit) = 
exp (—iHst/fh) p(t) exp (iHgt/h), we obtain the exact re- 
duced density operator in the matrix form of 



Pit) 



Poo(0) Poi(0)e- F (*) 

Pio(0)e- F *0 Pll (0) 



(18) 



with F(t) = icost + J drD{r). The same result was ob- 
tained in Ref. [19J | using the stochastic Schrodinger equa- 
tion approach. 

With the exact time evolution of the reduced density 
matrix, the one-time expectation value of the system op- 
erators 

(A (ti)) - Tr 5e i? [A (h) p T (0)] = Tr s [A (0) p (ti)] , 

(19) 
can be calculated exactly, where A(t\) represents a 
general system Heisenberg operator (s) and p{t\) — 
Trfl[/9j'(ti)] is the reduced Schrodinger density matrix 
operator at time t\ . We may also write in the interaction 
picture, 



{A(t 1 ))=Tr S(SR \A(t 1 )p T (t 1 ) 



= Tr, 



A(ti)p(h) 

(20) 



where pr is defined in Eq. {7J, p(t) = Tr R [pT {t)] and 
A(t) = exp(iH t/ft)Aexp(-iH t/K), and A = A(0). 

For a general system operator A = I , , 1 , we obtain 

exactly from either Eq. (JTSJ) or Eq. (J2U]) 

(A (h)) = e-fo 1 drD(r *> (ap w e l " stl + bpoie^^) 

+c P oo(0) + d Pll (0). (21) 

B. Two-time correlation functions 

In contrast to the Markovian case in which the QRT 
is valid, the time evolution of the reduced density matrix 



of a non-Markovian open system alone is not sufficient 
to obtain the two-time system operator CF's. This can 
be understood as follows. The two-time CF's of system 
operators A(ti)B(t 2 ) for t\ > t 2 can be written as 

(A{ti)B(t 2 )) 
= Tr s@R [U^(t 1 ,0)AU(t 1 ,0)U^t2,0)BU(t2, 0)p T (0)] 

= Trs®B.[AU(t 1 ,t 2 )BU(t 2l 0)p T (0)U\t 2l 0)U\t 1 ,tm) 

where the Heisenberg evolution operators U(ti,t 2 ) = 
U(h,0)W(t 2 ,0) and U(t,0) = exp(-iHt/h). If the 
environment is Markovian so the environment opera- 
tor CF at two different times is 5 correlated in time, 
then we may regard that the environment operator in 
U(ti,t 2 ) is not correlated with that in U(t 2 ,0). So the 
trace over the environment degrees of freedom for op- 
erator U(ti 7 t 2 ) and operator U(t 2 ,0) can be performed 
independently or separately. Thus one may first trace 
Pr(t 2 ) = U(t 2 ,0)pT{0)W(t 2 ,0) over the environment de- 
grees of freedom to obtain the reduced density matrix 
p(t 2 ) = TiR[pT(t 2 )}- Equation (|22l) in this case can be 
written as 

(A(ti)B(ta)) = Trs®R[AU(h,t2)(Bp(t 2 ) ® Ro)^ (t l7 t 2 )} 
= Tr s [A X (t)1 (23) 

where x( r ) is the effective reduced density matrix at time 
t = t\ — t 2 with the initial condition x(0) = Bp{t 2 ). 
Thus knowing the time evolution of the reduced density 
matrix in the Markovian case, one is able to calculate the 
two-time CF's of the system operators. This is also the 
reason why the QRT works in the Markovian case. But 
the situation differs for a non-Markovian environment as 
the environment operator in U(ti,t 2 ) may, in general, be 
correlated with that in t/(i 2 ,0). 

The two-time CF's of the system operators for the 
pure-dephasing spin-boson model can also be evaluated 
exactly. To evaluate the two-time CF of system opera- 
tors A(ti)B(t 2 ) for t\ > t 2 , we express it in terms of the 
interaction picture operators as 

(A(h)B(t 2 )) 
= Tr SBR lUHt 1 )AU(t 1 )U\t 2 )BU(t 2 )p T (0)} 
= TT S(BR [uHt 1 )A(t 1 )U(t 1 )U\t 2 )B(t 2 )U(t 2 )p T (0)l 

(24) 

where again an operator with a tilde on the top indi- 
cates that it is an operator in the interaction picture 
with respect to the free Hamiltonian Ho. Compared with 
Eq. (|2TH) . Eq. (|2"4"|) for general non-Markovian open sys- 
tems can not be expressed as a product of the reduced 
density matrix and system operators. So again, the re- 
duced density matrix alone is not sufficient to obtain the 
non-Markovian two-time system operator CF's. 

As we want to compare the results by the direct evalu- 
ation with those by the evolution equation ([3]), we calcu- 
late, in the following, the two-time CF's (A(ti)B{t 2 )) for 
different cases of system operators A and B. The struc- 
ture of the evolution equations in Ref. [15[ or Eqs. © and 



Q in this article depends on the commutation relations 
of operator A and operator L (or I/), and on the commu- 
tation relations of operator B and operator L(t — t 2 ) (or 
U{t — t 2 )), where L(t) = exp(iHst/H) Lexp(—iHst/h) 
is the system operator in the interaction picture with re- 
spect to Hs- For the pure-dephasing spin-boson model, 
Hs = (hu>s/2)a z , L — a z — L\ and then L'(t) = <r z . 
So we will discuss the two-time CF's in the follow- 
ing three cases and the trivial case of (a z (ti)a z (t 2 )} = 
(cr z (0)cr z (0)) = I is obvious due to [a z ,H] = 0. 

Case 1. [A,L] ^ and [B,L(t)] = 0. In this 
case, let us set A = aa + + ba~. and B — a z . Then 



A(t) 



aa+e 



iust 



bo-e- %ust and B(t) = a z . It is 



easy to see from Eq. (fT0)) that U(t) commutes with B{t) 
but anticommutes with A(t), i.e., [U(t),B(t)} = and 
{U(t),A(t)} — 0. Using these results and the fact that 
£/tM(i) [/M(£) = J ) we obtain from Eq. (|M|) 

(Ati)B(t 2 )) 

= -ap lt) (Oy^Tr R {U^ Ht 1 )m i Ht 1 )R Q } 

+bp ol (0)e^^Tv R [U^ 1 Ht 1 )U^(t 1 )R Q }. (25) 

Substituting the result of Eq. (TT5|) into Eq. (l25l , we arrive 
at the exact two-time CF's 



(Mh)B{t 2 )) 

p -J^drD(r) , 



-api O (0)e 



iusti 



bp i(0)e 



-iusti^ 



(26) 



Case 2. [A, L] = and [B, L(t)} ^ 0. In this case, let 
A = a z , and B = aa + + &<r_ . Similar to the calculations 
in Case 1, we obtain 



(A(h)B{t 2 )) 

e-f? dTD{T) (ap w (0)e^ - 6p O i(0)e 



— icJS*2^ 



(27) 

The exact two-time CF's of Eqs. (gB]) and (J2ZJ) depend on 
only one time variable, t\ or t 2 , respectively, since one of 
the system operator cr z (t) = a z (0) is time- independent. 

Case 3. [A,L] ^ and [B,L(t)] ^ 0. Suppose A = 
aa + + ba-, and B = a'a + + b'a-. In this case, both 
A{t) and B{t) anticommute with both U{t) and U^{t). 
Furthermore, A{t{)B{t 2 ) — ab'<r + <j~ exp[«ws(ii — ^2)] + 
6a'(T_cr + exp[— itusiti — t 2 )\. Thus we can obtain from 
Eq. dH 

(A(h)B(t 2 )) 
= a6'p 00 (0)e^ s(tl - t2) 

xTr^tm {h) (j{i} (t x )WW (t 2 )UW ( t2 )R a ] 

+toVii(0)e"" s(fl ~ f2) 
xTiRlU^Ht^U^Ht^U^Ht^U^Ht^Ro}. (28) 



It is obvious from Eq. (1281) that to evaluate the general 
two-time CF, we need to take into account the correla- 
tions of the reservoir operators of the evolution operators 



between different time periods of [0,£2] and [0,ii] be- 
fore the trace over the environment is performed. Using 
Eqs. (H2|), dHJ), and dHJ), we get 



Tt«[C^ t0} (*i)^ tl} (*i)CTt {1} ( fe ) J7* 0} (* a )i2o] 



exp 



where 



(ItD(t) 



<ItD(t) 



dTD(T,t 2 ) 



ft 2 

Jo 



D(r,t 2 ) — 4 / dsa c s(r - s). 



(29) 



(30) 



The term L 1 dT-D(r, £2) in Eq. (|2U)l describes the cross- 
time contribution of the environment CF's of the reser- 
voir operators in the evolution operators U^ n '(ti) and 
U Hn} (t 2 ) [or W^(h) and l> {n} (£ 2 )] of the two differ- 
ent time periods [0,ii] and [O,^]- We can see this from 
D(r,h) of Eq. ([30]) and in Eq. ([29]) that the environment 
CF a ff (t — s), defined in Eq. (IT7|) . has the time variable 
r in [0, ti] and the time variable s in [0, £2]- O n the other 
hand, the time evolution of the reduced density matrix 
fliTl) is involved with the reservoir operator CF's in the 
evolution operators of only the same time interval. As 
a result, it, alone, cannot provide us with the full infor- 
mation to evaluate the non-Markovian two-time CF, even 
in the weak system-environment coupling case. Similarly, 
we find that Tr R [UJ {1} (h)UJ 0} (h)^ (t 2 )UJ l} (t 2 )R Q ] 
has the same result as Eq. ([2"9"1) . Substituting these re- 
sults into Eq. (|28|) . finally we arrive at the two-time CF 



{A{h)B(t 2 )) 

r r l1 

exp — / cItD(t) 



dTD(r)+ J drD{TM) 



L Jo Jo Jo 

x (ab'p 00 (0)e^ s{tl - t2) +ba' Pll (0)e-^ s{tl - t2) 



(31) 



This non-Markovian finite-temperature two-time CF, to 
our knowledge, has not been presented in the literature. 



IV. EVALUATION BY DERIVED 

NON-MARKOVIAN FINITE-TEMPERATURE 

EVOLUTION EQUATIONS 

In this section, we will use the derived evolution equa- 
tions in Ref. [151 ] to compute the one-time expectation 
values and two-time CF's to compare with the exact ex- 
pressions evaluated in Sec. IIIII Despite the fact that the 
evolution equations in Ref. [15| derived perturbatively, 
the results obtained this way for the pure-dephasing spin- 
boson model happen to be the same as the exact expres- 
sions by the direct evaluation. 



A. Quantum master equation and one-time 
expectation values 

Before going to calculate the CF's, it is instructive to 
derive the master equation of the reduced system den- 
sity matrix for the model. After some calculations, we 
obtain for the Hamiltonian in the form of Eq. (QT) a time- 
covolutionless non-Markovian master equation [5|, [a, [2J- 
|28J ] valid to second order in the system-environment in- 
teraction strength 

dr{a(t - t)[L^L(t - t)p{t) - L(t - £)p(i)I f ] 

+(3{t - t)[LL\t - t)p{t) -D(t- t)p{t)L] 
+H.c.}, (32) 

where a{t — r) and (3{t — t) are defined in Eqs. @ and 
(O respectively, H.c. indicates the Hermitian conjugate 
of previous terms, and an operator with a tilde on the top 
indicates that it is an operator in the interaction picture. 
For the pure-dephasing spin-boson model, Eq. (|3"2"|) gives 
the master equation of the reduced system density matrix 



dp(t) 
dt 



-lUS 



?z,P(t)] 



D(t) 



[p{t)-a z p(t)a z ], (33) 



where D(t) is defined in Eq. (TT61) . It is not difficult 
to show that the exact expression of the density ma- 
trix (fT5f is the solution of the master equation (|3"31 al- 
though the master equation is derived perturbatively. 
Non-Markovian dynamics usually means that the cur- 
rent time evolution of the system state depends on its 
history, and the memory effects typically enters through 
integrals over the past state history. However, the non- 
Markovian system dynamics of some class of open quan- 
tum system models may be summed up and expressed 
as a time-local, convolutionless form [29] where the dy- 
namics is determined by the system state at the current 
time t only. This time-local, convolutionless class of open 
quantum systems may be treated exactly without any ap- 
proximation. The quantum Brownian motion model or 
the damped harmonic oscillator bilinearly coupled to a 
bosonic bath of harmonic oscillators |29l - [3l| is a famous 
example of this class. The pure-dephasing spin-boson 
model considered here also belongs to this class, and the 
non-Markovian effect in the master equation ([33]) is taken 
into account by the time-dependent coefficient D(t) in- 
stead of memory integral. This time-local, convolution- 
less property and the fact of [L, H s ] = allow the ex- 
act system density matrix Eq. (| 18|) to be obtained from 

Eq. U32D. 

Since the exact solution of the system density matrix 
([TH]) can be calculated from the perturbatively derived 
master equation (|33p . one may expect that the exact 
non-Markovian finite-temperature one-time expectation 
values and two-time CF's of the pure-dephasing model 
can be obtained from the evolution equation Q We show 



below that this is indeed the case, and at the same time 
the agreement of the results demonstrates the validity 
and practical usage of the evolution equation (j3|). 

For the pure-dephasing spin-boson model, H$ = 
(fiujs/2)a z , L = a z = L\ and we have Lft'(t) = a z . 
Taking A = <ji, % = x,y, z, in Eq. ([2]), we obtain straight- 
forwardly the evolution equations of the single-time ex- 
pectation values as 

d(a x (t 1 ))/dt 1 = -D(t 1 )(<r x (t 1 ))-uJs(cry(h)), (34) 
d (<r v (h)) Mi = -D(h) (o-„(ti)> + los M*i)) , (35) 
d(a z {ti))/dh = (36) 

with D{ti) defined in Eq. <JT6j ■ With proper chosen val- 
ues for a, b, c, and d of a general operator A for (Tj, one 
can verify that the exact expression of the expectation 
value of <Ji{ti) in Eq. flST) satisfies Eqs. ([34"j) -(136l). 



B. Two-time correlation functions 

Before using Eq. ([3]) to calculate the non-Markovian 
finite-temperature two-time CF's, we discuss briefly be- 
low the relation between the QRT and the evolution 
equation ([3]). If the last two terms of Eq. © vanish, 
then the single-time and two-time evolution equations 
@ and ([3J will have the same form with the same evo- 
lution coefficients and thus the QRT will be applicable. 
The last two terms of Eq. J3j or more generally the last 
term of Eq. (17) in Ref. [15J involve(s) the propagation 
from t = to r = t 2 , and these terms would vanish 
for the CF's (A(ti)B(O)) as t 2 = in this case. So the 
QRT is valid to calculate the CF's (A(t)B(0)) of both 
Markovian and non-Markovian open systems, where the 
system-environment density matrix is separable at t = 0. 
The QRT is also valid and is often applied to calcu- 
late, in the Markovian weak system-environment cou- 
pling case, more general CF's (A(t2 + r)B(t 2 )) or equiva- 
lent^ (A(ii)B(i 2 )) with t 2 j= 0. For example, the QRT is 
used to calculate the Markovian steady-state CF's, and 
in this case t 2 is set to any of the large times when the 
steady state is reached. This is because in the Marko- 
vian case, the last two terms of Eq. (J3j vanish since the 
time integration of the corresponding ^-correlated reser- 
voir CF's, a{t\ — t) oc <5(ti — r) and (3(ti — r) ex S(ti — r), 
over the variable r in the domain from to t 2 is zero as 
t\ > t 2 . On the other hand, the QRT cannot be blindly 
applied to calculate (A(ti)B(t 2 )) with t 2 ^ in a general 
non-Markovian open system due to the non-vanishing 
contributions of the cross correlation of the reservoir op- 
erators at two different times: a later time t\ and an 
earlier time in the period between and t 2 (see the last 
two terms of Eq. (J3J and also Fig. [1]). In other words, 
in contrast to the Markovian case, not only the initial 
condition (A(t 2 )B(t 2 )) for the two-time evolution equa- 
tion ([3]) but also the equation (J3J itself may depend on 
the choice of the starting time t 2 of the non-Markovian 
finite-temperature two-time CF's. In the steady state, 



the situation may change when t 2 is in any of the large 
times where the state and system expectation values do 
not change with time any more. In this case, the contri- 
butions from the last two terms of Eq. ([3J) saturate and 
do not depend on where time t 2 is set in the steady state, 
and thus both the Markovian and non-Markovian CF's 
may depend only on the time difference (t± —t 2 ) (see also 
Fig. [U . But the nonvanishing contributions from the last 
two terms of Eq. (J3j) would still make the non-Markovian 
CF's deviate from that obtained wrongly using the QRT 
in the non-Markovian case or obtained using the QRT in 
the Markovian case (see also Fig. [3J). 

For the time evolutions of system two-time CF's of the 
pure-dephasing spin-boson model, we also consider the 
following three cases as in Sec. IIIII Note that L(t) — 
a z = L\t). 

Case 1. [A,L] ^ and [B,L(i)] — 0. In this case, let 
A = ai,i = x, y and B = a z . By using Eq. ([3j), it is easy 
to obtain 



d(a x (ti)a z (t 2 ))/dti 

d(a y (t 1 )a z (t 2 ))/dti 



-D(ti)(o- B (*i)o-,(t 2 )> 
-u s {<r y (ti)<r z {t2)), (37) 

-D(tlX0j,(tlK(t2)) 
+V 8 ((T x {tl)(T t (t2)). (38) 



In this case, one can see that the evolution equations 
of single-time expectation values (<ji(ti)) t Eqs. (|M|) and 
(|35|). have the same forms as the evolution equations of 
two-time CF's (cri( i i) cr z( i 2)), Eqs. ([37J and ([gg|). respec- 
tively. Hence the QRT is valid in this case. It is easy to 
check that taking the derivative of Eq. ([25)1 with respect 
to t\ with a = b = 1 ( i.e. A= a x ) and a = —b = —i 
( i.e. A= a y ), one can obtain the evolution equations 
for (<r x {ti)a z (t 2 )) and (cr y {ti)a z (t 2 )), exactly the same 
as Eqs. (1371) and (135)) . respectively. 

Case 2. [A, L] = and [B, L(t)] ^ 0. In this case, let 
A = <j z and B — ai,i — x,y. By using Eq. ([3J, we then 
easily obtain 



d(<r z (ti)a x (t 2 ))/dti =0, 
d(a z (ti)a v (t 2 ))/dh =0. 



(39) 
(40) 

Indeed, Eq. fl^TJ) satisfies Eqs. (J3JJ and (gDJ), and 
{(J z (ti)a t (t 2 )) = {a z (t 2 )a t (t 2 )), independent of t%. 

Case 3. [A, L{t)] ^ and [B, L{t)] ^ 0. In this case, 
let A = ai, i = x, y and B — <jj,j — x,y. Eq. © straight- 
forwardly yields 

d(a x (ti)ay(t 2 )}/dh = -I>(ti)<a- B (*i)o- B (t a )> 

-u s (a y (tx)a y (t 2 )) -D(ti,t 2 )(a v (ti)a x (t 2 )), (41) 
d{a y {t l )a x {t 2 ))/dt l = -D(ti)(a y (tx)a x (t 2 )) 

+us(<rm(ti)a x (t2)) -D(h,t 2 )(a x (h)a y (h)), (42) 
d(a x (h)a x (t 2 )}/dh = -£>(ti)<a- x (*i)cr B (*2)> 

-u s (a y (ti)a x (t 2 )) + D(ti,t 2 )(a y (ti)a y (t 2 )), (43) 
d(a y (ti)a y (t 2 ))/dti = -D(ti)(a y (h)a y (t 2 )) 

+oj s {<T x {t 1 )oy{t 2 )) + £>(ii,i2)M*iV*(*2)), (44) 



where D(t\,ti) is defined in Eq. (l30l) . The evolution 
equations, Eqs. (|H ]) -(|13 |) . have different forms as those 
of single-time expectation values due to the existence of 
D(ti,t2) terms. As a result, the QRT does not hold in 
this case. Again, taking the derivative of Eq. (l3TT) with 
respect to t\ with properly chosen values for a, b, a' and 
&', we arrive at the same evolution equations as those 
from Eqs. (|4"Tj) to (|4"4"1) . Alternatively, solving the coupled 
equations, Eqs. (l4T1) - (|44p . one would obtain the solutions 
in a form as Eq. pip . 

The agreement between the results obtained by the di- 
rect operator evaluation and those obtained by solving 
the coupled evolution equations demonstrates clearly the 
validity of the equations © and ([3]) , In addition, the eas- 
iness to obtain Eqs. (j3~4")) -(f3l> | from the evolution equa- 
tion (0), and to obtain Eqs. (|37 ]l -([38 )l . Eqs. <(39j) and (|to|. 
and Eqs. (f4*T j) - (|4"4")) from the evolution equation §5§ illus- 
trates the practical usage of the non-Markovian finite- 
temperature evolution equations © and ([3]). 



V. RESULTS AND DISCUSSIONS 

To calculate the two important functions D(t) and 
D{t\,t 2 ), we need to evaluate the environment CF 



a c ff(£i 



duiJ(ui) {coth (hui/2kBT) cos[ui(ti — r)] 



-i siii[io(ti — t)]} 



(45) 



where J(u>) — J^ A \g\\ 2 S(uj — u>\) is the spectral density of 
the environment. We may consider any spectral density 
to characterize the environment, but for simplicity we 
consider an ohmic bath with exponential cut-off function 

as 



J{uj) = 7wexp(— w/A), 



(46) 



where A is the cut-off frequency and 7 is a dimensionless 
constant characterizing the interaction strength to the 
environment. At the zero temperature, the function D(t) 
and D(t\,t 2 ) have simple analytical forms: 



D(ti) = 4 7 



A 2 ^ 
1 + KH\ ' 



(47) 



4 7 A 2 t 2 [1 - KH.it, - t 2 ) - iA(2t! - h)] 
D{tlM) ~ (1 + A^Ml + A^-^] ^ 8) 

Consequently, the one-time expectation values and the 
two-time CF's also have simple analytical expressions. 
For example, the zero-temperature two-time CF's of 
Eq. (|26j) in case 1 and Eq. (|3"Tj) in case 3 are 



(A(h)B(t 2 )) 

(1 + A 2 t 2 )- 2 ^ (-a PlQ (0)e^ + bp 01 (0)e 
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FIG. 1. (Color online) Time evolutions of the real part 
of the system operator CF (a x (ti)a y (t2)) for four differ- 
ent cases: Markovian (solid line), non-Markovian using the 
QRT (dashed line) and non-Markovian (dot-dashed line) us- 
ing Eq. (J3]) and exact operator evaluation (dotted line). Other 
parameters used are cjs = 1, (JzbT /K) — 0.1, A = 5, 7 = 0.1, 
and £2 = 0.2. 



and 

(Mh)B(t 2 )) 

= [l + A 2 (h-t 2 ) 2 }- 2 -< 

4-yil arctan(A(ii— i2))+arctan(Ai2) %s J 

x (a6'poo(0)e l(JJs(tl ' t2) +6aVn(0)e"^ s(tl " t2) ) ,(50) 

respectively. 

Figure [T] shows the the time evolutions of the real part 
of the system operator CF (a x (t i)cT y (t 2 )) at a finite tem- 
perature of (hsT/h) = Q.Iujs- The time evolutions of the 
CF {<J x {ti)<J y {t 2 )) are obtained in four different cases: 
the first is in the Markovian case, the second is in the 
non-Markovian case with a finite cut-off frequency but 
wrongly applying the QRT method [i.e., neglecting the 
last two terms of Eq. (|3|) or equivalently neglecting the 
terms with D(ti,fa) in Eqs. (|II ]) -([ii |) ]. the third is in 
the non-Markovian case using the derived evolution equa- 
tions (pnj) - (T4"4"l) . and the fourth is the exact model using 
the direct operator evaluation. The initial environment 
state is in the thermal state and the system state is arbi- 
trarily chosen to be |*} = (^- |e) + | \g) J. The Marko- 
vian case in Fig. [Tj is described as follows. With the 
finite cut-off environment spectral density and with the 
system parameters used in Fig. [TJ the Markovian approx- 
imation may actually not be valid. If, however, we still 
assume that the environment correlation time in Eq. (|45j) 
is much smaller than all of the system time scales (i.e., 
Markovian approximation) , then we may replace the up- 
per time integration limit in Eq. (|16[) to infinity (i.e., 
t — > 00). As a result, the Markovian master equation 
or evolution equations can be obtained by just replac- 
ing the time-dependent coefficient D{t) in Eq. (|33p or 



8 °- 2 

& 
-0.1 



"^5-^ 






1,=0 

- -V - 1 

...« 2 =0.2 

t 2 =0.5 

..t 2 =l 

< 2 = 2 
.....« 2 =3 

t,-4 

■-V* 

---t 2 =10 


— Exact 

\ 






2 4 6 8 10 


l V — 






2 4 6 8 10 




— Markovian QRT 
- - - non-Mar kovian QRT 
— non-Markovian 
— Exact 



^rt 



t=t r l 2 



FIG. 2. (Color online) Time evolutions of the real part of the 
system operator CF {(T x (ti)<T y (t2)) for different values of £2 
using non-Markovian Eq. ((3]). The results of the time evolu- 
tions coincide with those obtained by the exact operator eval- 
uation. Other parameters used are cjs = 1, {ksT/K) — 0.1, 
A = 5, 7 = 0.1. The insets show the time evolutions of the 
real part of the expectation values (a x (t)) and {a y (t)). 



FIG. 3. (Color online) Time evolutions of the real part 
of the system operator CF {a x (ti)a y (t2)) for four differ- 
ent cases: Markovian (solid line), non-Markovian using the 
QRT (dashed line) and non-Markovian (dot-dashed line) us- 
ing Eq. (J3]) and exact operator evaluation (dotted line). Other 
parameters used are cjs = 1, (feT/S) = 0.1, A = 5, 7 = 0.1, 
and £2 = 10. 



in Eqs. (|34[) and ([55]) by its long-time limit value. At 
a finite temperature, the Markovian (time-independent) 
coefficient from Eqs. (|16l) . (fl5j) and P5|) can be written 
as 



£>oo = lim D{t) = 4~/irk B T/h. 

t—^-OO 



(51) 



We may see that Drx, — > as the temperature T — > 0. 
This is because in the Markovian limit, the decoherence 
or dephasing is strongly dependent on the infrared be- 
havior (w — > modes) of the environment in the pure 
dephasing model. Since the spectral density considered 
in Eq. (pjf5|) is Ohmic, we then have J(ui —> 0) = 0, and 
thus Deo — > at T = 0. This is in contrast to other 
quantum open system models with a resonant type of 
system-environment coupling, in which the environment 
modes near the system resonance frequency are relevant 
to the relaxation and decoherence. We can see from Fig.Q] 
that the difference between the results of the Markovian 
QRT case and the non-Markovian QRT case is visible, 
while the two-time CF's obtained by the non-Markovian 
evolution equation Q and by the exact operator eval- 
uation are identical for the pure-dephasing spin-boson 
model. The perfect agreement of the results between the 
non-Markovian evolution equation case and the exact op- 
erator evaluation case, and the significant difference in 
the short time region between the non-Markovian evolu- 
tion equation case and the wrong non-Markovian QRT 
case demonstrate clearly the validity and practical us- 
age of the evolution equation ([3]). All of the four cases 
approach one another to zero in the long time region. 

Figure [2] investigates the dependence of the exact two- 
time system operator CF on the time variable £2. We 
see that the time evolutions of the real part of the CF 
(c a; (ii)a'j / (i2)) as a function of t = £1 — t<i for the values of 



£2 < 1 behave quite differently, but they approach one an- 
other for £2 > 2. When £2 > 5, the steady state is reached 
as indicated in the time evolutions of the expectation val- 
ues (a x (t)) and (cr y (t)) shown in the insets of Fig. [2] In 
this case, the time evolutions of the two-time CF are 
independent of the choices of the starting time of £2 in 
the steady state and depend only on the time difference 
£ = £1 — £2 for the parameters used in Fig. [2] This can also 
be seen from the analytical expression of the exact zero- 
temperature CF (|50|) . The zero-temperature CF (J50| is 
a function of variables £1 and £2, but for a large value of 
the cut-off frequency A, when £2 is reasonably large, the 
CF depends almost only on £ = t\ — £2. Figure [3] shows 
the time evolutions of the real part of the steady-state 
(£2 = 10) system operator CF {<T x (ti)a y (t2)) obtained in 
four different cases as in Fig. [TJ As expected, the non- 
Markovian evolution equation case coincides with the ex- 
act operator evaluation case. They are, however, signifi- 
cantly different from the wrong non-Markovian QRT case 
and Markovian QRT case, even though the time evolu- 
tions of the steady-state two-time CF's depend only on 
the time difference £ = t\ — £2. One can also see that 
the CF's of the Markovian QRT and the non-Markovian 
QRT cases approach each other much more closely in the 
steady state than in Fig. [TJ 



VI. CONCLUSION 

We have evaluated the exact non-Markovian finite- 
temperature one-time expectation values and two-time 
CF's of the system operators for the exactly solvable 
pure-dephasing spin-boson model. The evaluation has 
been performed in two ways, one by exact direct opera- 
tor technique without any approximation and the other 
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by the evolution equations ([2]) and (|3]) valid to second or- 
der in the system-environment interaction Hamiltonian. 
Since the non-Markovian dynamics of the pure-dephasing 
spin-boson model can be cast into a time-local, convo- 
lutionless form and [L, H s ) — 0, the results obtained by 
the second-order evolution equations ([2]) and ([3]) turn out 
to be exactly the same as the exact results obtained by 
the exact direct operator evaluation. The agreement of 
the results between the two different approaches demon- 
strates clearly the validity of the evolution equations @ 
and ©. Furthermore, it is easy to obtain Eqs. (J3"T)F ([3"5 |) . 
Eqs. (gni) and ([30]). and Eqs. dHJ-dHJ) from the evolu- 
tion equation ([3]). Other non-Markovian open quantum 
system models that are not exactly solvable can be pro- 
ceeded in a similar way to obtain the time evolutions of 
their two-time system operator CF's valid to second or- 
der in the system-environment interaction Hamiltonian. 
This illustrates the practical usage of the evolution equa- 
tions. It is thus believed that the evolution equations @ 
and ([3]), which generalize the QRT to the non-Markovian 
finite-temperature case will have applications in many 
different branches of physics. 
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Appendix A: Derivation of time evolution operator 

To show the time evolution operator of Eq. (fTU|) . we 
begin from Eq. (|5|) with Hj(t) given by Eq. ([§]). Since 
Hi(t) in Eq. ^ contains only two major terms, which 



are, respectively, proportional to a\ and a\, one is tempt- 
ing to evaluate the time-ordered exponent by the reverse 
operator identity of Eq. (fl"2"|) 



(Al) 



e A+B = e A e B e'^ A ' B \ 



valid for the commutator [A, B] commuting with both A 
and B. As the exponent operates at different times and 
[Hj(t),Hi(r)] ^ 0, it is not correct to use the precise 
form of Eq. (|A1[) . The proper procedure done in [20, |32j 
is to separate the two terms of Eq. © in the time-ordered 
exponent of Eq. (J5J by 




i J* dr^gxvze^aU 



d Te i lo ds £ a 9x <r z e"°* s a{ 



^9\<J z e- luJxT a x ) e-'^^^w^^M 



(A2) 

where T is the time-ordering operator. Then using the 
identity e~^ a ^a\e^ a >- — a\ + <\> in the exponent of the 
time-ordered term in Eq. (|A2[) . we obtain 



U(t) = exp 



•p(-*/ dT^ 5A ^e-^ r a A ) 



t l-T 



exp - / dr dsJ2 \9\\ 2 e~ tUx{T ~ s) A3) 



We have dropped the time-ordering operator in Eq. (|A3|) . 
Using the operator identity of Eq. (fT2|) to combine 
the first two terms in Eq. (|A3|) . we then obtain 
Eq. (flO"]). Note that the first term in Eq. flTUJ is just 

exp[— (i/%) J drHilr)] if the time-ordering operator for 
U(t) is not performed. The correct time-ordering proce- 
dure generates extra phase factors in Eq. (|A3|) and thus 
inEq. (fTTJ|) . 
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